home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / DBio / ureadseq.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-17  |  59.4 KB  |  2,173 lines  |  [TEXT/R*ch]
  1. /* File: ureadseq.c
  2.  *
  3.  * Reads and writes nucleic/protein sequence in various
  4.  * formats. Data files may have multiple sequences.
  5.  *
  6.  * Copyright 1990 by d.g.gilbert
  7.  * biology dept., indiana university, bloomington, in 47405
  8.  * e-mail: gilbertd@bio.indiana.edu
  9.  *
  10.  * This program may be freely copied and used by anyone.
  11.  * Developers are encourged to incorporate parts in their
  12.  * programs, rather than devise their own private sequence
  13.  * format.
  14.  *
  15.  * This should compile and run with any ANSI C compiler.
  16.  *
  17.  */
  18.  
  19.  
  20. #include <stdio.h>
  21. #include <ctype.h>
  22. #include <string.h>
  23. #include <stdlib.h>
  24.  
  25. #define UREADSEQ_G
  26. #include "ureadseq.h"
  27.  
  28.  
  29.  
  30. int Strcasecmp(const char *a, const char *b)  /* from Nlm_StrICmp */
  31. {
  32.   int diff, done;
  33.   if (a == b)  return 0;
  34.   done = 0;
  35.   while (! done) {
  36.     diff = to_upper(*a) - to_upper(*b);
  37.     if (diff) return diff;
  38.     if (*a == '\0') done = 1;
  39.     else { a++; b++; }
  40.     }
  41.   return 0;
  42. }
  43.  
  44. int Strncasecmp(const char *a, const char *b, long maxn) /* from Nlm_StrNICmp */
  45. {
  46.   int diff, done;
  47.   if (a == b)  return 0;
  48.   done = 0;
  49.   while (! done) {
  50.     diff = to_upper(*a) - to_upper(*b);
  51.     if (diff) return diff;
  52.     if (*a == '\0') done = 1;
  53.     else {
  54.       a++; b++; maxn--;
  55.       if (! maxn) done = 1;
  56.       }
  57.     }
  58.   return 0;
  59. }
  60.  
  61.  
  62.  
  63. #    define MALLOC(size)            (char*)malloc(size)
  64. #    define REALLOC(p,size)     p= (char*)realloc(p,size)
  65.  
  66. #ifndef Local
  67. # define Local      static    /* local functions */
  68. #endif
  69.  
  70. #define kStartLength 500     /* 500000 to temp fix Unix bug */
  71.  
  72. const char *aminos      = "ABCDEFGHIKLMNPQRSTVWXYZ*";
  73. const char *primenuc    = "ACGTU";
  74. const char *protonly    = "EFIPQZ";
  75.  
  76. const char kNocountsymbols[5]  = "_.-?";
  77. const char stdsymbols[6]  = "_.-*?";
  78. const char allsymbols[32] = "_.-*?<>{}[]()!@#$%^&=+;:'/|`~\"\\";
  79. static const char *seqsymbols   = allsymbols;
  80.  
  81. const char nummask[11]   = "0123456789";
  82. const char nonummask[11] = "~!@#$%^&*(";
  83.  
  84. FILE    * gFile;
  85. long    gLinestart = 0;
  86.  
  87.  
  88.  
  89.  
  90. /*
  91.     use general form of isseqchar -- all chars + symbols.
  92.     no formats except nbrf (?) use symbols in data area as
  93.     anything other than sequence chars.
  94. */
  95.  
  96.  
  97.  
  98.                           /* Local variables for readSeq: */
  99. struct ReadSeqVars {
  100.   short choice, err, nseq;
  101.   long  seqlen, maxseq, seqlencount;
  102.   short topnseq;
  103.   long  topseqlen;
  104.   const char *fname;
  105.   char *seq, *seqid, matchchar;
  106.   boolean allDone, done, filestart, addit;
  107.   FILE  *f;
  108.   long  linestart;
  109.   char  s[256], *sp;
  110.  
  111.   int (*isseqchar)(int);  
  112.     ReadWriteProc    reader;
  113. };
  114.  
  115.  
  116.  
  117. #define EOFreader(V)    (boolean)(*V->reader)((char *)NULL, 0L, kRSFile_Eof)
  118. #define REWINDreader(V)    (*V->reader)((char *)NULL, 0L, kRSFile_Rewind)
  119.  
  120. Local long  readWriteFromFile(char* line, long maxline, short action)
  121. {
  122.     static boolean geof = false;
  123.     long   bytes;
  124.     /*fprintf(stderr," r/w action= %d,  maxline= %ld ", action, maxline);*/
  125.     /* ! one case where _Read returns 0 bytes at end of file, but _Eof returns false !! */
  126.     
  127.     switch (action) {
  128.     
  129.     case kRSFile_Eof:
  130.         if (geof) {  /* fix for bad feof() ! */
  131.             char ch= fgetc( gFile);
  132.             if (ch != EOF) { ungetc( ch, gFile); geof= false; }
  133.             }
  134.         else geof= feof( gFile);
  135.         return geof;
  136.  
  137.     case kRSFile_Read:
  138.         bytes= (long) fgets (line, (int)maxline, gFile);
  139.         if (!bytes) geof= true; /* fix for bad feof() ! */
  140.         return bytes;
  141.     
  142.     case kRSFile_Write:
  143.         geof= false;
  144.         fputs( line, gFile);
  145.         return 0;
  146.  
  147.     case kRSFile_Seek:
  148.         geof= false;
  149.         fseek( gFile, maxline, 0);
  150.         return 0;
  151.  
  152.     case kRSFile_Rewind:
  153.         rewind(gFile);
  154.         geof= false;
  155.         return 0;
  156.  
  157.     case kRSFile_Tell:
  158.         return ftell(gFile);
  159.         
  160.     default:
  161.         return 0;
  162.     }
  163. }
  164.  
  165.  
  166.  
  167.  
  168. int isSeqChar(int c)
  169. {
  170.  /*  return (isalpha(c) || strchr(seqsymbols,c)); */
  171.   if (isalpha(c) || strchr(seqsymbols,c)) return c;
  172.   else return 0;
  173. }
  174.  
  175. int isGCGSeqChar(int c)
  176. {
  177.  /*  return (isalpha(c) || strchr(seqsymbols,c)); */
  178.   if (isalpha(c) || strchr(seqsymbols,c))  {
  179.     if (c == '.') return '-'; /* do the indel translate */
  180.     else return c;
  181.     }
  182.   else return 0;
  183. }
  184.  
  185. int isSeqNumChar(int c)
  186. {
  187.   /* return (isalnum(c) || strchr(seqsymbols,c)); */
  188.   if (isalnum(c) || strchr(seqsymbols,c)) return c;
  189.   else return 0;
  190. }
  191.  
  192.  
  193. int isAnyChar(int c)
  194. {
  195.  /*  return isprint(c); /* wrap in case isascii is macro */
  196.  if (isprint(c)) return c;
  197.  else return 0;
  198. }
  199.  
  200.  
  201. Local void callreadline( ReadWriteProc reader, char *s)
  202. {
  203.   char  *cp;
  204.  
  205.     /*fprintf(stderr,"callreadline &reader= %ld \n", reader);*/
  206.     
  207.     gLinestart= (*reader)( (char*)NULL, 0L, kRSFile_Tell);
  208.     if (0 == (*reader)( s, 256L, kRSFile_Read) )
  209.         *s = 0;
  210.   else {
  211.     cp = strchr(s, '\n');
  212.     if (cp != NULL) *cp = 0;
  213.     }
  214. }
  215.  
  216. /********
  217. Local void readline(FILE *f, char *s, long *linestart)
  218. {
  219.   char  *cp;
  220.  
  221.     *linestart= ftell(f);
  222.     if (NULL == fgets(s, 256, f))
  223.     *s = 0;
  224.   else {
  225.     cp = strchr(s, '\n');
  226.     if (cp != NULL) *cp = 0;
  227.     }
  228. }
  229. ******/
  230.  
  231. Local void getline(struct ReadSeqVars *V)
  232. {
  233.   /*readline(V->f, V->s, &V->linestart);  */
  234.     callreadline(V->reader, V->s);
  235. }
  236.  
  237. Local void ungetline(struct ReadSeqVars *V)
  238. {
  239.     /*fseek(V->f, V->linestart, 0); */
  240.     (void) (*(V->reader))( (char*)NULL, gLinestart, kRSFile_Seek);  
  241. }
  242.  
  243.  
  244. Local void addseq(char *s, struct ReadSeqVars *V)
  245. {
  246.   char  *ptr;
  247.   register char seqc;
  248.   
  249.   if (V->addit) while (*s != 0) {
  250.     if ((seqc= (V->isseqchar)(*s)) != 0) {
  251.       if (V->seqlen >= V->maxseq) {
  252.         V->maxseq += kStartLength;
  253. #ifdef NOTmacintosh
  254.                 SetPtrSize( (Ptr)V->seq, V->maxseq+1);
  255.                 if (MemError() != 0) {
  256.           V->err = eMemFull;
  257.           return;
  258.           }
  259. #else
  260.         ptr = (char*) realloc(V->seq, V->maxseq+1);
  261.         if (ptr==NULL) {
  262.           V->err = eMemFull;
  263.           return;
  264.           }
  265.         else V->seq = ptr;
  266. #endif
  267.                 }
  268.       V->seq[(V->seqlen)++] = seqc; /* *s */
  269.       }
  270.     s++;
  271.     }
  272. }
  273.  
  274. Local void countseq(char *s, struct ReadSeqVars *V)
  275.  /* this must count all valid seq chars, for some formats (paup-sequential) even
  276.     if we are skipping seq... */
  277. {
  278.   while (*s != 0) {
  279.     if ((V->isseqchar)(*s)) {
  280.       (V->seqlencount)++;
  281.       }
  282.     s++;
  283.     }
  284. }
  285.  
  286.  
  287. Local void addinfo(char *s, struct ReadSeqVars *V)
  288. {
  289.   char s2[256], *si;
  290.   boolean saveadd;
  291.   int (*oldIsSeqChar)(int);  
  292.   
  293.   si = s2;
  294.   while (*s == ' ') s++;
  295.   sprintf(si, " %d)  %s\n", V->nseq, s);
  296.  
  297.   saveadd = V->addit;
  298.   V->addit = true;
  299.   oldIsSeqChar= V->isseqchar;
  300.   V->isseqchar = isAnyChar;
  301.   addseq( si, V);
  302.   V->addit = saveadd;
  303.   V->isseqchar = oldIsSeqChar;  
  304. }
  305.  
  306.  
  307.  
  308.  
  309. Local void readLoop(short margin, boolean addfirst,
  310.             boolean (*endTest)(boolean *addend, boolean *ungetend, struct ReadSeqVars *V),
  311.             struct ReadSeqVars *V)
  312. {
  313.   boolean addend = false;
  314.   boolean ungetend = false;
  315.  
  316.   V->nseq++;
  317.   if (V->choice == kListSequences) V->addit = false;
  318.   else V->addit = (V->nseq == V->choice);
  319.   if (V->addit) V->seqlen = 0;
  320.  
  321.   if (addfirst) addseq(V->s, V);
  322.   do {
  323.     getline(V);
  324.     V->done = EOFreader(V); /* feof(V->f); */
  325.     V->done |= (*endTest)( &addend, &ungetend, V);
  326.     if (V->addit && (addend || !V->done) && (strlen(V->s) > margin)) {
  327.       addseq( (V->s)+margin, V);
  328.     }
  329.   } while (!V->done);
  330.  
  331.   if (V->choice == kListSequences) addinfo(V->seqid, V);
  332.   else {
  333.     V->allDone = (V->nseq >= V->choice);
  334.     if (V->allDone && ungetend) ungetline(V);
  335.     }
  336. }
  337.  
  338.  
  339.  
  340. Local boolean endIG( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  341. {
  342.   *addend = true; /* 1 or 2 occur in line w/ bases */
  343.   *ungetend= false;
  344.   return((strchr(V->s,'1')!=NULL) || (strchr(V->s,'2')!=NULL));
  345. }
  346.  
  347. Local void readIG(struct ReadSeqVars *V)
  348. {
  349. /* 18Aug92: new IG format -- ^L between sequences in place of ";" */
  350.   char  *si;
  351.  
  352.   while (!V->allDone) {
  353.     do {
  354.       getline(V);
  355.       for (si= V->s; *si != 0 && *si < ' '; si++) *si= ' '; /* drop controls */
  356.       if (*si == 0) *V->s= 0; /* chop line to empty */
  357.     } while (! (EOFreader(V) || ((*V->s != 0) && (*V->s != ';') ) ));
  358.     if (EOFreader(V))
  359.       V->allDone = true;
  360.     else {
  361.       strcpy(V->seqid, V->s);
  362.       readLoop(0, false, endIG, V);
  363.       }
  364.   }
  365. }
  366.  
  367.  
  368.  
  369. Local boolean endStrider( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  370. {
  371.   *addend = false;
  372.   *ungetend= false;
  373.   return (strstr( V->s, "//") != NULL);
  374. }
  375.  
  376. Local void readStrider(struct ReadSeqVars *V)
  377. { /* ? only 1 seq/file ? */
  378.  
  379.   while (!V->allDone) {
  380.     getline(V);
  381.     if (strstr(V->s,"; DNA sequence  ") == V->s)
  382.       strcpy(V->seqid, (V->s)+16);
  383.     else
  384.       strcpy(V->seqid, (V->s)+1);
  385.     while ((!EOFreader(V)) && (*V->s == ';')) {
  386.       getline(V);
  387.       }
  388.     if (EOFreader(V)) V->allDone = true;
  389.     else readLoop(0, true, endStrider, V);
  390.   }
  391. }
  392.  
  393.  
  394. Local boolean endPIR( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  395. {
  396.   *addend = false;
  397.   *ungetend= (strstr(V->s,"ENTRY") == V->s);
  398.   return ((strstr(V->s,"///") != NULL) || *ungetend);
  399. }
  400.  
  401. Local void readPIR(struct ReadSeqVars *V)
  402. { /*PIR -- many seqs/file */
  403.  
  404.   while (!V->allDone) {
  405.     while (! (EOFreader(V) 
  406.                 || strstr(V->s,"ENTRY")  || strstr(V->s,"SEQUENCE")) )
  407.       getline(V);
  408.     strcpy(V->seqid, (V->s)+16);
  409.     while (! (EOFreader(V) || strstr(V->s,"SEQUENCE") == V->s))
  410.       getline(V);
  411.     readLoop(0, false, endPIR, V);
  412.  
  413.     if (!V->allDone) {
  414.      while (! (EOFreader(V) || ((*V->s != 0)
  415.        && (strstr( V->s,"ENTRY") == V->s))))
  416.         getline(V);
  417.       }
  418.     if (EOFreader(V)) V->allDone = true;
  419.   }
  420. }
  421.  
  422.  
  423. Local boolean endGB( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  424. {
  425.   *addend = false;
  426.   *ungetend= (strstr(V->s,"LOCUS") == V->s);
  427.   return ((strstr(V->s,"//") != NULL) || *ungetend);
  428. }
  429.  
  430. Local void readGenBank(struct ReadSeqVars *V)
  431. { /*GenBank -- many seqs/file */
  432.  
  433.   while (!V->allDone) {
  434.     strcpy(V->seqid, (V->s)+12);
  435.     while (! (EOFreader(V) || strstr(V->s,"ORIGIN") == V->s))
  436.       getline(V);
  437.     readLoop(0, false, endGB, V);
  438.  
  439.     if (!V->allDone) {
  440.      while (! (EOFreader(V) || ((*V->s != 0)
  441.        && (strstr( V->s,"LOCUS") == V->s))))
  442.         getline(V);
  443.       }
  444.     if (EOFreader(V)) V->allDone = true;
  445.   }
  446. }
  447.  
  448.  
  449. Local boolean endNBRF( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  450. {
  451.   char  *a;
  452.  
  453.   if ((a = strchr(V->s, '*')) != NULL) { /* end of 1st seq */
  454.     /* "*" can be valid base symbol, drop it here */
  455.     *a = 0;
  456.     *addend = true;
  457.     *ungetend= false;
  458.     return(true);
  459.     }
  460.   else if (*V->s == '>') { /* start of next seq */
  461.     *addend = false;
  462.     *ungetend= true;
  463.     return(true);
  464.     }
  465.   else
  466.     return(false);
  467. }
  468.  
  469. Local void readNBRF(struct ReadSeqVars *V)
  470. {
  471.   while (!V->allDone) {
  472.     strcpy(V->seqid, (V->s)+4);
  473.     getline(V);   /*skip title-junk line*/
  474.     readLoop(0, false, endNBRF, V);
  475.     if (!V->allDone) {
  476.      while (!(EOFreader(V) || (*V->s != 0 && *V->s == '>')))
  477.         getline(V);
  478.       }
  479.     if (EOFreader(V)) V->allDone = true;
  480.   }
  481. }
  482.  
  483.  
  484.  
  485. Local boolean endPearson( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  486. {
  487.   *addend = false;
  488.   *ungetend= true;
  489.   return(*V->s == '>');
  490. }
  491.  
  492. Local void readPearson(struct ReadSeqVars *V)
  493. {
  494.   while (!V->allDone) {
  495.     strcpy(V->seqid, (V->s)+1);
  496.     readLoop(0, false, endPearson, V);
  497.     if (!V->allDone) {
  498.      while (!(EOFreader(V) || ((*V->s != 0) && (*V->s == '>'))))
  499.         getline(V);
  500.       }
  501.     if (EOFreader(V)) V->allDone = true;
  502.   }
  503. }
  504.  
  505.  
  506.  
  507. Local boolean endEMBL( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  508. {
  509.   *addend = false;
  510.   *ungetend= (strstr(V->s,"ID   ") == V->s);
  511.   return ((strstr(V->s,"//") != NULL) || *ungetend);
  512. }
  513.  
  514. Local void readEMBL(struct ReadSeqVars *V)
  515. {
  516.   while (!V->allDone) {
  517.     strcpy(V->seqid, (V->s)+5);
  518.     do {
  519.       getline(V);
  520.     } while (!(EOFreader(V) || (strstr(V->s,"SQ   ") == V->s)));
  521.  
  522.     readLoop(0, false, endEMBL, V);
  523.     if (!V->allDone) {
  524.       while (!(EOFreader(V) ||
  525.          ((*V->s != '\0') && (strstr(V->s,"ID   ") == V->s))))
  526.       getline(V);
  527.     }
  528.     if (EOFreader(V)) V->allDone = true;
  529.   }
  530. }
  531.  
  532.  
  533.  
  534. Local boolean endZuker( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  535. {
  536.   *addend = false;
  537.   *ungetend= true;
  538.   return( *V->s == '(' );
  539. }
  540.  
  541. Local void readZuker(struct ReadSeqVars *V)
  542. {
  543.   /*! 1st string is Zuker's Fortran format */
  544.  
  545.   while (!V->allDone) {
  546.     getline(V);  /*s == "seqLen seqid string..."*/
  547.     strcpy(V->seqid, (V->s)+6);
  548.     readLoop(0, false, endZuker, V);
  549.     if (!V->allDone) {
  550.       while (!(EOFreader(V) |
  551.         ((*V->s != '\0') && (*V->s == '('))))
  552.           getline(V);
  553.       }
  554.     if (EOFreader(V)) V->allDone = true;
  555.   }
  556. }
  557.  
  558.  
  559.  
  560. Local boolean endFitch( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  561. {
  562.   /* this is a somewhat shaky end,
  563.     1st char of line is non-blank for seq. title
  564.   */
  565.   *addend = false;
  566.   *ungetend= true;
  567.   return( *V->s != ' ' );
  568. }
  569.  
  570. Local void readFitch(struct ReadSeqVars *V)
  571. {
  572.   boolean first;
  573.  
  574.   first = true;
  575.   while (!V->allDone) {
  576.     if (!first) strcpy(V->seqid, V->s);
  577.     readLoop(0, first, endFitch, V);
  578.     if (EOFreader(V)) V->allDone = true;
  579.     first = false;
  580.     }
  581. }
  582.  
  583.  
  584. Local void readPlain(struct ReadSeqVars *V)
  585. {
  586.   V->nseq++;
  587.   V->addit = (V->choice > 0);
  588.   if (V->addit) V->seqlen = 0;
  589.   addseq(V->seqid, V);   /*from above..*/
  590.   if (V->fname!=NULL) sprintf(V->seqid, "%s  [Unknown form]", V->fname);
  591.   else sprintf(V->seqid, "  [Unknown form]");
  592.   do {
  593.     addseq(V->s, V);
  594.     V->done = EOFreader(V);
  595.     getline(V);
  596.   } while (!V->done);
  597.   if (V->choice == kListSequences) addinfo(V->seqid, V);
  598.   V->allDone = true;
  599. }
  600.  
  601.  
  602. Local void readUWGCG(struct ReadSeqVars *V)
  603. {
  604. /*
  605. 10nov91: Reading GCG files casued duplication of last line when
  606.          EOF followed that line !!!
  607.     fix: getline now sets *V->s = 0
  608. */
  609.   char  *si;
  610.  
  611.   V->nseq++;
  612.   V->addit = (V->choice > 0);
  613.   if (V->addit) V->seqlen = 0;
  614.   strcpy(V->seqid, V->s);
  615.   /*writeseq: "    %s  Length: %d  (today)  Check: %d  ..\n" */
  616.   /*drop above or ".." from id*/
  617.   if ((si = strstr(V->seqid,"  Length: ")) != 0) *si = 0;
  618.   else if ((si = strstr(V->seqid,"..")) != 0) *si = 0;
  619.   do {
  620.     V->done = EOFreader(V);
  621.     getline(V);
  622.     if (!V->done) addseq((V->s), V);
  623.   } while (!V->done);
  624.   if (V->choice == kListSequences) addinfo(V->seqid, V);
  625.   V->allDone = true;
  626. }
  627.  
  628.  
  629. Local void readOlsen(struct ReadSeqVars *V)
  630. { /* G. Olsen /print output from multiple sequence editor */
  631.  
  632.   char    *si, *sj, *sk, *sm, sid[40], snum[20];
  633.   boolean indata = false;
  634.   short snumlen = 0;
  635.  
  636.   V->addit = (V->choice > 0);
  637.   if (V->addit) V->seqlen = 0;
  638.   REWINDreader(V); 
  639.     V->nseq= 0;
  640.   do {
  641.     getline(V);
  642.     V->done = EOFreader(V);
  643.  
  644.     if (V->done && !(*V->s)) break;
  645.     else if (indata) {
  646.       if ( (si= strstr(V->s, sid))
  647.         && ( (sm= strstr(V->s, snum)) != NULL ) && (sm < si - snumlen) ) {
  648.         /* && (strstr(V->s, snum) == si - snumlen - 1) ) */
  649.  
  650.         /* Spaces are valid alignment data !! */
  651. /* 17Oct91: Error, the left margin is 21 not 22! */
  652. /* dropped some nucs up to now -- my example file was right shifted ! */
  653. /* variable right id margin, drop id-2 spaces at end */
  654. /*
  655.   VMS CC COMPILER (VAXC031) mess up:
  656.   -- Index of 21 is chopping 1st nuc on VMS systems Only!
  657.   Byte-for-byte same ame rnasep.olsen sequence file !
  658. */
  659.  
  660.         /* si = (V->s)+21; < was this before VMS CC wasted my time */
  661.         si += 10;  /* use strstr index plus offset to outfox VMS CC bug */
  662.  
  663.         if ((sk = strstr(si, sid))!=0) *(sk-2) = 0;
  664.         for (sk = si; *sk != 0; sk++) {
  665.            if (*sk == ' ') *sk = '.';
  666.            /* 18aug92: !! some olsen masks are NUMBERS !! which addseq eats */
  667.            else if (isdigit(*sk)) *sk= nonummask[*sk - '0'];
  668.            }
  669.  
  670.         addseq(si, V);
  671.         }
  672.       }
  673.  
  674.     else if ((sk = strstr(V->s, "): "))!=0) {  /* seq info header line */
  675.   /* 18aug92: correct for diff seqs w/ same name -- use number, e.g. */
  676.   /*   3 (Agr.tume):  agrobacterium.prna  18-JUN-1987 16:12 */
  677.   /* 328 (Agr.tume):  agrobacterium.prna XYZ  19-DEC-1992   */
  678.       (V->nseq)++;
  679.       si = 1 + strchr(V->s,'(');
  680.       *sk = ' ';
  681.       if (V->choice == kListSequences) addinfo( si, V);
  682.       else if (V->nseq == V->choice) {
  683.         strcpy(V->seqid, si);
  684.         sj = strchr(V->seqid, ':');
  685.         while (*(--sj) == ' ') ;
  686.         while (--sj != V->seqid) { if (*sj == ' ') *sj = '_'; }
  687.  
  688.         *sk = 0;
  689.         while (*(--sk) == ' ') *sk = 0;
  690.         strcpy(sid, si);
  691.  
  692.         si= V->s;
  693.         while ((*si <= ' ') && (*si != 0)) si++;
  694.         snumlen=0;
  695.         while (si[snumlen] > ' ' && snumlen<20)
  696.          { snum[snumlen]= si[snumlen]; snumlen++; }
  697.         snum[snumlen]= 0;
  698.         }
  699.  
  700.       }
  701.  
  702.     else if (strstr(V->s,"identity:   Data:")) {
  703.       indata = true;
  704.       if (V->choice == kListSequences) V->done = true;
  705.       }
  706.  
  707.   } while (!V->done);
  708.  
  709.   V->allDone = true;
  710. } /*readOlsen*/
  711.  
  712.  
  713. Local void readMSF(struct ReadSeqVars *V)
  714. { /* gcg's MSF, mult. sequence format, interleaved ! */
  715.  
  716.   char    *si, *sj, sid[128];
  717.   boolean indata = false;
  718.   int     atseq= 0, iline= 0;
  719.  
  720.   V->addit = (V->choice > 0);
  721.   if (V->addit) V->seqlen = 0;
  722.   REWINDreader(V); 
  723.     V->nseq= 0;
  724.   do {
  725.     getline(V);
  726.     V->done = EOFreader(V);
  727.  
  728.     if (V->done && !(*V->s)) break;
  729.     else if (indata) {
  730.       /*somename  ...gpvedai .......t.. aaigr..vad tvgtgptnse aipaltaaet */
  731.       /*       E  gvenae.kgv tentna.tad fvaqpvylpe .nqt...... kv.affynrs */
  732.  
  733.       si= V->s;
  734.       skipwhitespace(si);
  735.       /* for (sj= si; isalnum(*sj); sj++) ; bug -- delwiche uses "-", "_" and others in names*/
  736.       for (sj= si; *sj > ' '; sj++) ;
  737.       *sj= 0;
  738.       if ( *si ) {
  739.         if ( (0==strcmp(si, sid)) ) {
  740.           addseq(sj+1, V);
  741.           }
  742.         iline++;
  743.         }
  744.       }
  745.  
  746.     else if (NULL != (si = strstr(V->s, "Name: "))) {  /* seq info header line */
  747.       /* Name: somename      Len:   100  Check: 7009  Weight:  1.00 */
  748.  
  749.       (V->nseq)++;
  750.       si += 6;
  751.       if (V->choice == kListSequences) addinfo( si, V);
  752.       else if (V->nseq == V->choice) {
  753.         strcpy(V->seqid, si);
  754.         si = V->seqid;
  755.         skipwhitespace(si);
  756.         /* for (sj= si; isalnum(*sj); sj++) ; -- bug */
  757.         for (sj= si; *sj > ' '; sj++) ;
  758.         *sj= 0;
  759.         strcpy(sid, si);
  760.         }
  761.       }
  762.  
  763.     else if ( strstr(V->s,"//") /*== V->s*/ )  {
  764.       indata = true;
  765.       iline= 0;
  766.       if (V->choice == kListSequences) V->done = true;
  767.       }
  768.  
  769.   } while (!V->done);
  770.  
  771.  
  772.   V->allDone = true;
  773. } /*readMSF*/
  774.  
  775.  
  776.  
  777. Local void readPAUPinterleaved(struct ReadSeqVars *V)
  778. { /* PAUP mult. sequence format, interleaved or sequential! */
  779.  
  780.   char    *si, *sj, *send, sid[40], sid1[40], saveseq[255];
  781.   boolean first = true, indata = false, domatch;
  782.   int     atseq= 0, iline= 0, ifmc, saveseqlen=0;
  783.  
  784. #define fixmatchchar(s) { \
  785.   for (ifmc=0; ifmc<saveseqlen; ifmc++) \
  786.     if (s[ifmc] == V->matchchar) s[ifmc]= saveseq[ifmc]; }
  787.  
  788.   V->addit = (V->choice > 0);
  789.   V->seqlencount = 0;
  790.   if (V->addit) V->seqlen = 0;
  791.   /* rewind(V->f); V->nseq= 0;  << do in caller !*/
  792.   indata= true; /* call here after we find "matrix" */
  793.   domatch= (V->matchchar > 0);
  794.  
  795.   do {
  796.     getline(V);
  797.     V->done = EOFreader(V);
  798.  
  799.     if (V->done && !(*V->s)) break;
  800.     else if (indata) {
  801.       /* [         1                    1                    1         ]*/
  802.       /* human     aagcttcaccggcgcagtca ttctcataatcgcccacggR cttacatcct*/
  803.       /* chimp     ................a.t. .c.................a ..........*/
  804.       /* !! need to correct for V->matchchar */
  805.       si= V->s;
  806.       skipwhitespace(si);
  807.       if (strchr(si,';')) indata= false;
  808.  
  809.       if (isalnum(*si))  {
  810.         /* valid data line starts w/ a left-justified seq name in columns [0..8] */
  811.         if (first) {
  812.           (V->nseq)++;
  813.           if (V->nseq >= V->topnseq) first= false;
  814.           for (sj = si; isalnum(*sj); sj++) ;
  815.           send= sj;
  816.           skipwhitespace(sj);
  817.           if (V->choice == kListSequences) {
  818.             *send= 0;
  819.             addinfo( si, V);
  820.             }
  821.           else if (V->nseq == V->choice) {
  822.             if (domatch) {
  823.               if (V->nseq == 1) { strcpy( saveseq, sj); saveseqlen= strlen(saveseq); }
  824.               else fixmatchchar( sj);
  825.               }
  826.             addseq(sj, V);
  827.             *send= 0;
  828.             strcpy(V->seqid, si);
  829.             strcpy(sid, si);
  830.             if (V->nseq == 1) strcpy(sid1, sid);
  831.             }
  832.           }
  833.  
  834.         else if ( (strstr(si, sid) == si) ){
  835.           while (isalnum(*si)) si++;
  836.           skipwhitespace(si);
  837.           if (domatch) {
  838.             if (V->nseq == 1) { strcpy( saveseq, si); saveseqlen= strlen(saveseq); }
  839.             else fixmatchchar( si);
  840.             }
  841.           addseq(si, V);
  842.           }
  843.  
  844.         else if (domatch && (strstr(si, sid1) == si)) {
  845.           strcpy( saveseq, si);
  846.           saveseqlen= strlen(saveseq);
  847.           }
  848.  
  849.         iline++;
  850.         }
  851.       }
  852.  
  853.     else if ( strstr(V->s,"matrix") )  {
  854.       indata = true;
  855.       iline= 0;
  856.       if (V->choice == kListSequences) V->done = true;
  857.       }
  858.  
  859.   } while (!V->done);
  860.  
  861.   V->allDone = true;
  862. } /*readPAUPinterleaved*/
  863.  
  864.  
  865.  
  866. Local void readPAUPsequential(struct ReadSeqVars *V)
  867. { /* PAUP mult. sequence format, interleaved or sequential! */
  868.   char    *si, *sj;
  869.   boolean atname = true, indata = false;
  870.  
  871.   V->addit = (V->choice > 0);
  872.   if (V->addit) V->seqlen = 0;
  873.   V->seqlencount = 0;
  874.   /* REWINDreader(V); V->nseq= 0;  << do in caller !*/
  875.   indata= true; /* call here after we find "matrix" */
  876.   do {
  877.     getline(V);
  878.     V->done = EOFreader(V);
  879.  
  880.     if (V->done && !(*V->s)) break;
  881.     else if (indata) {
  882.       /* [         1                    1                    1         ]*/
  883.       /* human     aagcttcaccggcgcagtca ttctcataatcgcccacggR cttacatcct*/
  884.       /*           aagcttcaccggcgcagtca ttctcataatcgcccacggR cttacatcct*/
  885.       /* chimp     ................a.t. .c.................a ..........*/
  886.       /*           ................a.t. .c.................a ..........*/
  887.  
  888.       si= V->s;
  889.       skipwhitespace(si);
  890.       if (strchr(si,';')) indata= false;
  891.       if (isalnum(*si))  {
  892.         /* valid data line starts w/ a left-justified seq name in columns [0..8] */
  893.         if (atname) {
  894.           (V->nseq)++;
  895.           V->seqlencount = 0;
  896.           atname= false;
  897.           sj= si+1;
  898.           while (isalnum(*sj)) sj++;
  899.           if (V->choice == kListSequences) {
  900.             /* !! we must count bases to know when topseqlen is reached ! */
  901.             countseq(sj, V);
  902.             if (V->seqlencount >= V->topseqlen) atname= true;
  903.             *sj= 0;
  904.             addinfo( si, V);
  905.             }
  906.           else if (V->nseq == V->choice) {
  907.             addseq(sj, V);
  908.             V->seqlencount= V->seqlen;
  909.             if (V->seqlencount >= V->topseqlen) atname= true;
  910.             *sj= 0;
  911.             strcpy(V->seqid, si);
  912.             }
  913.           else {
  914.             countseq(sj, V);
  915.             if (V->seqlencount >= V->topseqlen) atname= true;
  916.             }
  917.           }
  918.  
  919.         else if (V->nseq == V->choice) {
  920.           addseq(V->s, V);
  921.           V->seqlencount= V->seqlen;
  922.           if (V->seqlencount >= V->topseqlen) atname= true;
  923.           }
  924.         else {
  925.           countseq(V->s, V);
  926.           if (V->seqlencount >= V->topseqlen) atname= true;
  927.           }
  928.         }
  929.       }
  930.  
  931.     else if ( strstr(V->s,"matrix") )  {
  932.       indata = true;
  933.       atname= true;
  934.       if (V->choice == kListSequences) V->done = true;
  935.       }
  936.  
  937.   } while (!V->done);
  938.  
  939.   V->allDone = true;
  940. } /*readPAUPsequential*/
  941.  
  942.  
  943. Local void readPhylipInterleaved(struct ReadSeqVars *V)
  944. {
  945.   char    *si, *sj;
  946.   boolean first = true;
  947.   int     iline= 0;
  948.  
  949.   V->addit = (V->choice > 0);
  950.   if (V->addit) V->seqlen = 0;
  951.   V->seqlencount = 0;
  952.   /* sscanf( V->s, "%d%d", &V->topnseq, &V->topseqlen); << topnseq == 0 !!! bad scan !! */
  953.   si= V->s;
  954.   skipwhitespace(si);
  955.   V->topnseq= atoi(si);
  956.   while (isdigit(*si)) si++;
  957.   skipwhitespace(si);
  958.   V->topseqlen= atol(si);
  959.   /* fprintf(stderr,"Phylip-ileaf: topnseq=%d  topseqlen=%d\n",V->topnseq, V->topseqlen); */
  960.  
  961.   do {
  962.     getline(V);
  963.     V->done = EOFreader(V);
  964.  
  965.     if (V->done && !(*V->s)) break;
  966.     si= V->s;
  967.     skipwhitespace(si);
  968.     if (*si != 0) {
  969.  
  970.       if (first) {  /* collect seq names + seq, as fprintf(outf,"%-10s  ",seqname); */
  971.         (V->nseq)++;
  972.         if (V->nseq >= V->topnseq) first= false;
  973.         sj= V->s+10;  /* past name, start of data */
  974.         if (V->choice == kListSequences) {
  975.           *sj= 0;
  976.           addinfo( si, V);
  977.           }
  978.         else if (V->nseq == V->choice) {
  979.           addseq(sj, V);
  980.           *sj= 0;
  981.           strcpy(V->seqid, si);
  982.           }
  983.         }
  984.       else if ( iline % V->nseq == V->choice -1 ) {
  985.         addseq(si, V);
  986.         }
  987.       iline++;
  988.     }
  989.   } while (!V->done);
  990.  
  991.   V->allDone = true;
  992. } /*readPhylipInterleaved*/
  993.  
  994.  
  995.  
  996. Local boolean endPhylipSequential( boolean *addend, boolean *ungetend, struct ReadSeqVars *V)
  997. {
  998.   boolean done;
  999.     /* *addend = false; << need this true if EOF(file) on last dataline */
  1000.   *ungetend= false;
  1001.   countseq( V->s, V);
  1002.     done= (V->seqlencount >= V->topseqlen);
  1003.   *addend = !done;
  1004.     return done;
  1005. }
  1006.  
  1007. Local void readPhylipSequential(struct ReadSeqVars *V)
  1008. {
  1009.   short  i;
  1010.   char  *si;
  1011.   /* sscanf( V->s, "%d%d", &V->topnseq, &V->topseqlen); < ? bad sscan ? */
  1012.   si= V->s;
  1013.   skipwhitespace(si);
  1014.   V->topnseq= atoi(si);
  1015.   while (isdigit(*si)) si++;
  1016.   skipwhitespace(si);
  1017.   V->topseqlen= atol(si);
  1018.   getline(V);
  1019.   while (!V->allDone) {
  1020.     V->seqlencount= 0;
  1021.     strncpy(V->seqid, (V->s), 10);
  1022.     V->seqid[10]= 0;
  1023.     for (i=0; i<10 && V->s[i]; i++) V->s[i]= ' ';
  1024.     readLoop(0, true, endPhylipSequential, V);
  1025.     if (EOFreader(V)) V->allDone = true;
  1026.     }
  1027. }
  1028.  
  1029.  
  1030.  
  1031.  
  1032. Local void readSeqMain(
  1033.       struct ReadSeqVars *V,
  1034.       const long  skiplines,
  1035.       const short format)
  1036. {
  1037. #define tolowerstr(s) { long Itlwr, Ntlwr= strlen(s); \
  1038.   for (Itlwr=0; Itlwr<Ntlwr; Itlwr++) s[Itlwr]= to_lower(s[Itlwr]); }
  1039.  
  1040.   boolean gotuw;
  1041.   long l;
  1042.  
  1043.     /*fprintf(stderr,"readSeqMain &reader= %ld  &gFile= %ld\n", V->reader, gFile);*/
  1044.  
  1045.   V->linestart= 0;
  1046.   V->matchchar= 0;
  1047.   if (V->reader == NULL)
  1048.     V->err = eFileNotFound;
  1049.   else {
  1050.  
  1051.     for (l = skiplines; l > 0; l--) getline( V);
  1052.  
  1053.     do {
  1054.       getline( V);
  1055.       for (l= strlen(V->s); (l > 0) && (V->s[l] == ' '); l--) ;
  1056.     } while ((l == 0) && !EOFreader(V));
  1057.  
  1058.     if (EOFreader(V)) V->err = eNoData;
  1059.  
  1060.     else switch (format) {
  1061.       case kPlain : readPlain(V); break;
  1062.       case kIG    : readIG(V); break;
  1063.       case kStrider: readStrider(V); break;
  1064.       case kGenBank: readGenBank(V); break;
  1065.       case kPIR   : readPIR(V); break;
  1066.       case kNBRF  : readNBRF(V); break;
  1067.       case kPearson: readPearson(V); break;
  1068.       case kEMBL  : readEMBL(V); break;
  1069.       case kZuker : readZuker(V); break;
  1070.       case kOlsen : readOlsen(V); break;
  1071.       case kMSF   : 
  1072.           V->isseqchar = isGCGSeqChar;
  1073.                 readMSF(V); 
  1074.                 break;
  1075.  
  1076.       case kPAUP    : {
  1077.         boolean done= false;
  1078.         boolean interleaved= false;
  1079.         char  *cp;
  1080.         /* REWINDreader(V); V->nseq= 0; ?? assume it is at top ?? skiplines ... */
  1081.         while (!done) {
  1082.           getline( V);
  1083.           tolowerstr( V->s);
  1084.           if (strstr( V->s, "matrix")) done= true;
  1085.           if (strstr( V->s, "interleav")) interleaved= true;
  1086.           if (NULL != (cp=strstr( V->s, "ntax=")) )  V->topnseq= atoi(cp+5);
  1087.           if (NULL != (cp=strstr( V->s, "nchar=")) )  V->topseqlen= atoi(cp+6);
  1088.           if (NULL != (cp=strstr( V->s, "matchchar=")) )  {
  1089.             cp += 10;
  1090.             if (*cp=='\'') cp++;
  1091.             else if (*cp=='"') cp++;
  1092.             V->matchchar= *cp;
  1093.             }
  1094.           }
  1095.         if (interleaved) readPAUPinterleaved(V);
  1096.         else readPAUPsequential(V);
  1097.         }
  1098.         break;
  1099.  
  1100.       /* kPhylip: ! can't determine in middle of file which type it is...*/
  1101.       /* test for interleave or sequential and use Phylip4(ileave) or Phylip2 */
  1102.       case kPhylip2:
  1103.         readPhylipSequential(V);
  1104.         break;
  1105.       case kPhylip4: /* == kPhylip3 */
  1106.         readPhylipInterleaved(V);
  1107.         break;
  1108.  
  1109.       default:
  1110.         V->err = eUnknownFormat;
  1111.         break;
  1112.  
  1113.       case kFitch :
  1114.         strcpy(V->seqid, V->s); getline(V);
  1115.         readFitch(V);
  1116.         break;
  1117.  
  1118.       case kGCG:
  1119.         V->isseqchar = isGCGSeqChar;
  1120.         do {
  1121.           gotuw = (strstr(V->s,"..") != NULL);
  1122.           if (gotuw) readUWGCG(V);
  1123.           getline(V);
  1124.         } while (!(EOFreader(V) || V->allDone));
  1125.         break;
  1126.       }
  1127.     }
  1128.  
  1129.   V->filestart= false;
  1130.   V->seq[V->seqlen] = 0; /* stick a string terminator on it */
  1131. }
  1132.  
  1133.  
  1134. char *readSeqFp(
  1135.       const short whichEntry,  /* index to sequence in file */
  1136.       FILE  *fp,   /* pointer to open seq file */
  1137.       const long  skiplines,
  1138.       const short format,      /* sequence file format */
  1139.       long  *seqlen,     /* return seq size */
  1140.       short *nseq,       /* number of seqs in file, for listSeqs() */
  1141.       short *error,      /* return error */
  1142.       char  *seqid)      /* return seq name/info */
  1143. {
  1144.   struct ReadSeqVars V;
  1145.  
  1146.     gFile= fp; gLinestart = 0;
  1147.     V.reader = readWriteFromFile;
  1148.     
  1149.   if (format < kMinFormat || format > kMaxFormat) {
  1150.     *error = eUnknownFormat;
  1151.     *seqlen = 0;
  1152.     return NULL;
  1153.     }
  1154.  
  1155.   V.choice = whichEntry;
  1156.   V.fname  = NULL;  /* don't know */
  1157.   V.seq    = (char*) MALLOC(kStartLength+1);
  1158.   V.maxseq = kStartLength;
  1159.   V.seqlen = 0;
  1160.   V.seqid  = seqid;
  1161.  
  1162.   V.f = fp;
  1163.   V.filestart= (ftell( fp) == 0); 
  1164.   /* !! in sequential read, must remove current seq position from choice/whichEntry counter !! ... */
  1165.   if (V.filestart)  V.nseq = 0;
  1166.   else V.nseq= *nseq;  /* track where we are in file...*/
  1167.  
  1168.   *V.seqid = '\0';
  1169.   V.err = 0;
  1170.   V.nseq = 0;
  1171.   V.isseqchar = isSeqChar;
  1172.   if (V.choice == kListSequences) ; /* leave as is */
  1173.   else if (V.choice <= 0) V.choice = 1; /* default ?? */
  1174.   V.addit = (V.choice > 0);
  1175.   V.allDone = false;
  1176.  
  1177.   readSeqMain(&V, skiplines, format);
  1178.  
  1179.   *error = V.err;
  1180.   *seqlen = V.seqlen;
  1181.   *nseq = V.nseq;
  1182.   return V.seq;
  1183. }
  1184.  
  1185. char *readSeq(
  1186.       const short whichEntry,  /* index to sequence in file */
  1187.       const char  *filename,   /* file name */
  1188.       const long  skiplines,
  1189.       const short format,      /* sequence file format */
  1190.       long  *seqlen,     /* return seq size */
  1191.       short *nseq,       /* number of seqs in file, for listSeqs() */
  1192.       short *error,      /* return error */
  1193.       char  *seqid)      /* return seq name/info */
  1194. {
  1195.   struct ReadSeqVars V;
  1196.  
  1197.   if (format < kMinFormat || format > kMaxFormat) {
  1198.     *error = eUnknownFormat;
  1199.     *seqlen = 0;
  1200.     return NULL;
  1201.     }
  1202.  
  1203.   V.choice = whichEntry;
  1204.   V.fname  = filename;  /* don't need to copy string, just ptr to it */
  1205.   V.seq    = (char*) MALLOC( kStartLength+1);
  1206.   V.maxseq = kStartLength;
  1207.   V.seqlen = 0;
  1208.   V.seqid  = seqid;
  1209.  
  1210.   V.f = NULL;
  1211.   *V.seqid = '\0';
  1212.   V.err = 0;
  1213.   V.nseq = 0;
  1214.   V.isseqchar = isSeqChar;
  1215.   if (V.choice == kListSequences) ; /* leave as is */
  1216.   else if (V.choice <= 0) V.choice = 1; /* default ?? */
  1217.   V.addit = (V.choice > 0);
  1218.   V.allDone = false;
  1219.  
  1220.   V.f = fopen(V.fname, "r");
  1221.   V.filestart= true;
  1222.     gFile= V.f; gLinestart = 0;
  1223.     V.reader = readWriteFromFile;
  1224.  
  1225.   readSeqMain(&V, skiplines, format);
  1226.  
  1227.   if (V.f != NULL) fclose(V.f);
  1228.   *error = V.err;
  1229.   *seqlen = V.seqlen;
  1230.   *nseq = V.nseq;
  1231.   return V.seq;
  1232. }
  1233.  
  1234.  
  1235.  
  1236. char *readSeqCall(
  1237.       const short whichEntry,  /* index to sequence in file */
  1238.         ReadWriteProc reader,
  1239.       const long  skiplines,
  1240.       const short format,      /* sequence file format */
  1241.       long  *seqlen,     /* return seq size */
  1242.       short *nseq,       /* number of seqs in file, for listSeqs() */
  1243.       short *error,      /* return error */
  1244.       char  *seqid)      /* return seq name/info */
  1245. {
  1246.   struct ReadSeqVars V;
  1247.  
  1248.   if (format < kMinFormat || format > kMaxFormat) {
  1249.     *error = eUnknownFormat;
  1250.     *seqlen = 0;
  1251.     return NULL;
  1252.     }
  1253.  
  1254.   V.choice = whichEntry;
  1255.   V.fname  = NULL;  /* don't know */
  1256.   V.seq    = (char*) MALLOC( kStartLength+1);
  1257.   V.maxseq = kStartLength;
  1258.   V.seqlen = 0;
  1259.   V.seqid  = seqid;
  1260.   *V.seqid = '\0';
  1261.  
  1262.   V.f = NULL;
  1263.   V.err = 0;
  1264.   V.nseq = 0;
  1265.   V.isseqchar = isSeqChar;
  1266.   if (V.choice == kListSequences) ; /* leave as is */
  1267.   else if (V.choice <= 0) V.choice = 1; /* default ?? */
  1268.   V.addit = (V.choice > 0);
  1269.   V.allDone = false;
  1270.  
  1271.     gFile= NULL; gLinestart = 0;
  1272.     V.reader = reader;
  1273.      V.filestart= (0 == (*reader)((char*)NULL, 0L, kRSFile_Tell));
  1274.  
  1275.   readSeqMain(&V, skiplines, format);
  1276.  
  1277.   *error = V.err;
  1278.   *seqlen = V.seqlen;
  1279.   *nseq = V.nseq;
  1280.   return V.seq;
  1281. }
  1282.  
  1283.  
  1284.  
  1285. char *listSeqs(
  1286.       const char  *filename,   /* file name */
  1287.       const long skiplines,
  1288.       const short format,      /* sequence file format */
  1289.       short *nseq,       /* number of seqs in file, for listSeqs() */
  1290.       short *error)      /* return error */
  1291. {
  1292.   char  seqid[256];
  1293.   long  seqlen;
  1294.  
  1295.     /*fprintf(stderr,"listSeqs filename= %s  format= %d  nseq= %d\n", filename, format, *nseq);*/
  1296.   return readSeq( kListSequences, filename, skiplines, format,
  1297.                   &seqlen, nseq, error, seqid);
  1298. }
  1299.  
  1300. char *listSeqsCall(
  1301.       ReadWriteProc    reader,
  1302.       const long skiplines,
  1303.       const short format,      /* sequence file format */
  1304.       short *nseq,       /* number of seqs in file, for listSeqs() */
  1305.       short *error)      /* return error */
  1306. {
  1307.   char  seqid[256];
  1308.   long  seqlen;
  1309.  
  1310.   return readSeqCall( kListSequences, reader, skiplines, format,
  1311.                   &seqlen, nseq, error, seqid);
  1312. }
  1313.  
  1314.  
  1315.  
  1316.  
  1317. short seqFileFormat(    /* return sequence format number, see ureadseq.h */
  1318.     const char *filename,
  1319.     long  *skiplines,   /* return #lines to skip any junk like mail header */
  1320.     short *error)       /* return any error value or 0 */
  1321. {
  1322.   FILE      *fseq;
  1323.   short      format;
  1324.  
  1325.   fseq  = fopen(filename, "r");
  1326.   format= seqFileFormatFp( fseq, skiplines, error);
  1327.   if (fseq!=NULL) fclose(fseq);
  1328.   return format;
  1329. }
  1330.  
  1331. short seqFileFormatFp(
  1332.     FILE *fseq,
  1333.     long  *skiplines,   /* return #lines to skip any junk like mail header */
  1334.     short *error)       /* return any error value or 0 */
  1335. {
  1336.     gFile = fseq; gLinestart = 0;
  1337.     return seqFileFormatCall( readWriteFromFile, skiplines, error);
  1338. }
  1339.  
  1340.         
  1341.  
  1342. short seqFileFormatCall(
  1343.     ReadWriteProc reader,
  1344.     long  *skiplines,   /* return #lines to skip any junk like mail header */
  1345.     short *error)       /* return any error value or 0 */
  1346. {
  1347.   boolean   foundDNA= false, foundIG= false, foundStrider= false,
  1348.             foundGB= false, foundPIR= false, foundEMBL= false, foundNBRF= false,
  1349.             foundPearson= false, foundFitch= false, foundPhylip= false, foundZuker= false,
  1350.             gotolsen= false, gotpaup = false, gotasn1 = false, gotuw= false, gotMSF= false,
  1351.             isfitch= false,  isphylip= false, done= false;
  1352.   short     format= kUnknown;
  1353.   int       nlines= 0, k, splen= 0, otherlines= 0, aminolines= 0, dnalines= 0;
  1354.   char      sp[256];
  1355.   long      linestart=0;
  1356.   int         maxlines2check=500;
  1357.  
  1358. #define ReadOneLine(sp)   \
  1359.          { done |= (boolean)((*reader)((char*)NULL, 0L, kRSFile_Eof)); \
  1360.          callreadline( reader, sp);  \
  1361.     if (!done) { splen = strlen(sp); ++nlines; } }
  1362.  
  1363. /*  { done |= feof(fp); \ */
  1364. /*      readline( fp, sp, &linestart); \ */
  1365.  
  1366.     /*fprintf(stderr,"seqFileFormatCall &reader= %ld  &gFile= %ld\n", reader, gFile); */
  1367.  
  1368.   *skiplines = 0;
  1369.   *error = 0;
  1370.   if (reader == NULL) { *error = eFileNotFound;  return kNoformat; }
  1371.  
  1372.   while ( !done ) {
  1373.     ReadOneLine(sp);
  1374.  
  1375.     /* check for mailer head & skip past if found */
  1376.     if (nlines < 4 && !done) {
  1377.       if ((strstr(sp,"From ") == sp) || (strstr(sp,"Received:") == sp)) {
  1378.         do {
  1379.           /* skip all lines until find one blank line */
  1380.           ReadOneLine(sp);
  1381.           if (!done) for (k=0; (k<splen) && (sp[k]==' '); k++) ;
  1382.           } while ((!done) && (k < splen));
  1383.         *skiplines = nlines; /* !? do we want #lines or #bytes ?? */
  1384.         }
  1385.       }
  1386.  
  1387.     if (sp==NULL || *sp==0)
  1388.       ; /* nada */
  1389.  
  1390.     /* high probability identities: */
  1391.  
  1392.     else if ( strstr(sp,"MSF:") && strstr(sp,"Type:") && strstr(sp,"Check:") )
  1393.       gotMSF= true;
  1394.  
  1395.     else if ((strstr(sp,"..") != NULL) && (strstr(sp,"Check:") != NULL))
  1396.       gotuw= true;
  1397.  
  1398.     else if (strstr(sp,"identity:   Data:") != NULL)
  1399.       gotolsen= true;
  1400.  
  1401.     else if ( strstr(sp,"::=") &&
  1402.       (strstr(sp,"Bioseq") ||       /* Bioseq or Bioseq-set */
  1403.        strstr(sp,"Seq-entry") ||
  1404.        strstr(sp,"Seq-submit") ) )  /* can we read submit format? */
  1405.           gotasn1= true;
  1406.  
  1407.     else if ( strstr(sp,"#NEXUS") == sp )
  1408.       gotpaup= true;
  1409.  
  1410.     /* uncertain identities: */
  1411.  
  1412.     else if (*sp ==';') {
  1413.       if (strstr(sp,"Strider") !=NULL) foundStrider= true;
  1414.       else foundIG= true;
  1415.       }
  1416.  
  1417.     else if (strstr(sp,"LOCUS") == sp)
  1418.       foundGB= true;
  1419.     else if (strstr(sp,"ORIGIN") == sp)
  1420.       foundGB= true;
  1421.  
  1422.     else if (strstr(sp,"ENTRY   ") == sp)  /* ? also (strcmp(sp,"\\\\\\")==0) */
  1423.       foundPIR= true;
  1424.     else if (strstr(sp,"SEQUENCE") == sp)
  1425.       foundPIR= true;
  1426.  
  1427.     else if (*sp == '>') {
  1428.       if (sp[3] == ';') foundNBRF= true;
  1429.       else foundPearson= true;
  1430.       }
  1431.  
  1432.     else if (strstr(sp,"ID   ") == sp)
  1433.       foundEMBL= true;
  1434.     else if (strstr(sp,"SQ   ") == sp)
  1435.       foundEMBL= true;
  1436.  
  1437.     else if (*sp == '(')
  1438.       foundZuker= true;
  1439.  
  1440.     else {
  1441.       if (nlines - *skiplines == 1) {
  1442.         int ispp= 0, ilen= 0;
  1443.         sscanf( sp, "%d%d", &ispp, &ilen);
  1444.         if (ispp > 0 && ilen > 0) isphylip= true;
  1445.         }
  1446.       else if (isphylip && nlines - *skiplines == 2) {
  1447.         int  tseq;
  1448.         tseq= getseqtype(sp+10, strlen(sp+10));
  1449.         if ( isalpha(*sp)     /* 1st letter in 2nd line must be of a name */
  1450.          && (tseq != kOtherSeq))  /* sequence section must be okay */
  1451.             foundPhylip= true;
  1452.         }
  1453.  
  1454.       for (k=0, isfitch= true; isfitch && (k < splen); k++) {
  1455.         if (k % 4 == 0) isfitch &= (sp[k] == ' ');
  1456.         else isfitch &= (sp[k] != ' ');
  1457.         }
  1458.       if (isfitch && (splen > 20)) foundFitch= true;
  1459.  
  1460.       /* kRNA && kDNA are fairly certain...*/
  1461.       switch (getseqtype( sp, splen)) {
  1462.         case kOtherSeq: otherlines++; break;
  1463.         case kAmino   : if (splen>20) aminolines++; break;
  1464.         case kDNA     :
  1465.         case kRNA     : if (splen>20) dnalines++; break;
  1466.         case kNucleic : break; /* not much info ? */
  1467.         }
  1468.  
  1469.       }
  1470.  
  1471.                     /* pretty certain */
  1472.     if (gotolsen) {
  1473.       format= kOlsen;
  1474.       done= true;
  1475.       }
  1476.     else if (gotMSF) {
  1477.       format= kMSF;
  1478.       done= true;
  1479.       }
  1480.     else if (gotasn1) {
  1481.       /* !! we need to look further and return  kASNseqentry | kASNseqset */
  1482.       /*
  1483.         seqentry key is Seq-entry ::=
  1484.         seqset key is Bioseq-set ::=
  1485.         ?? can't read these yet w/ ncbi tools ??
  1486.           Seq-submit ::=
  1487.           Bioseq ::=  << fails both bioseq-seq and seq-entry parsers !
  1488.       */
  1489.       if (strstr(sp,"Bioseq-set")) format= kASNseqset;
  1490.       else if (strstr(sp,"Seq-entry")) format= kASNseqentry;
  1491.       else format= kASN1;  /* other form, we can't yet read... */
  1492.       done= true;
  1493.       }
  1494.     else if (gotpaup) {
  1495.       format= kPAUP;
  1496.       done= true;
  1497.       }
  1498.  
  1499.     else if (gotuw) {
  1500.       if (foundIG) format= kIG;  /* a TOIG file from GCG for certain */
  1501.       else format= kGCG;
  1502.       done= true;
  1503.       }
  1504.  
  1505.     else if ((dnalines > 1) || done || (nlines > maxlines2check)) {
  1506.           /* decide on most likely format */
  1507.           /* multichar idents: */
  1508.       if (foundStrider) format= kStrider;
  1509.       else if (foundGB) format= kGenBank;
  1510.       else if (foundPIR) format= kPIR;
  1511.       else if (foundEMBL) format= kEMBL;
  1512.       else if (foundNBRF) format= kNBRF;
  1513.           /* single char idents: */
  1514.       else if (foundIG) format= kIG;
  1515.       else if (foundPearson) format= kPearson;
  1516.       else if (foundZuker) format= kZuker;
  1517.           /* digit ident: */
  1518.       else if (foundPhylip) format= kPhylip;
  1519.           /* spacing ident: */
  1520.       else if (foundFitch) format= kFitch;
  1521.           /* no format chars: */
  1522.       else if (otherlines > 0) format= kUnknown;
  1523.       else if (dnalines > 1) format= kPlain;
  1524.       else if (aminolines > 1) format= kPlain;
  1525.       else format= kUnknown;
  1526.  
  1527.       done= true;
  1528.       }
  1529.  
  1530.     /* need this for possible long header in olsen format */
  1531.      else if (strstr(sp,"): ") != NULL)
  1532.        maxlines2check++;
  1533.     }
  1534.  
  1535.   if (format == kPhylip) {
  1536.     /* check for interleaved or sequential -- really messy */
  1537.     int tname, tseq;
  1538.     long i, j, nspp= 0, nlen= 0, ilen, leaf= 0, seq= 0;
  1539.     char  *ps;
  1540.  
  1541.         (void) (*reader)((char*)NULL, 0L, kRSFile_Rewind);
  1542.         /*rewind(fp);*/
  1543.     for (i=0; i < *skiplines; i++) ReadOneLine(sp);
  1544.     nlines= 0;
  1545.     ReadOneLine(sp);
  1546.     sscanf( sp, "%d%d", &nspp, &nlen);
  1547.     ReadOneLine(sp); /* 1st seq line */
  1548.     for (ps= sp+10, ilen=0; *ps!=0; ps++) if (isprint(*ps)) ilen++;
  1549.  
  1550.     for (i= 1; i<nspp; i++) {
  1551.       ReadOneLine(sp);
  1552.  
  1553.       tseq= getseqtype(sp+10, strlen(sp+10));
  1554.       tname= getseqtype(sp, 10);
  1555.       for (j=0, ps= sp; isspace(*ps) && j<10; ps++, j++)  ;
  1556.       for (ps= sp; *ps!=0; ps++) if (isprint(*ps)) ilen++;
  1557.  
  1558.       /* find probable interleaf or sequential ... */
  1559.       if (j>=9) seq += 10; /* pretty certain not ileaf */
  1560.       else {
  1561.         if (tseq != tname) leaf++; else seq++;
  1562.         if (tname == kDNA || tname == kRNA) seq++; else leaf++;
  1563.         }
  1564.  
  1565.       if (ilen <= nlen && j<9) {
  1566.         if (tname == kOtherSeq) leaf += 10;
  1567.         else if (tname == kAmino || tname == kDNA || tname == kRNA) seq++; else leaf++;
  1568.         }
  1569.       else if (ilen > nlen) {
  1570.         ilen= 0;
  1571.         }
  1572.       }
  1573.     for ( nspp *= 2 ; i<nspp; i++) {  /* this should be only bases if interleaf */
  1574.       ReadOneLine(sp);
  1575.  
  1576.       tseq= getseqtype(sp+10, strlen(sp+10));
  1577.       tname= getseqtype(sp, 10);
  1578.       for (ps= sp; *ps!=0; ps++) if (isprint(*ps)) ilen++;
  1579.       for (j=0, ps= sp; isspace(*ps) && j<10; ps++, j++)  ;
  1580.       if (j<9) {
  1581.         if (tname == kOtherSeq) seq += 10;
  1582.         if (tseq != tname) seq++; else leaf++;
  1583.         if (tname == kDNA || tname == kRNA) leaf++; else seq++;
  1584.         }
  1585.       if (ilen > nlen) {
  1586.         if (j>9) leaf += 10; /* must be a name here for sequent */
  1587.         else if (tname == kOtherSeq) seq += 10;
  1588.         ilen= 0;
  1589.         }
  1590.       }
  1591.  
  1592.     if (leaf > seq) format= kPhylip4;
  1593.     else format= kPhylip2;
  1594.     }
  1595.  
  1596.     /*fprintf(stderr, "seqFileFormatCall &reader= %ld  &gFile= %ld  format= %ld\n", reader, gFile, format);*/
  1597.   return(format);
  1598. #undef  ReadOneLine
  1599. } /* SeqFileFormat */
  1600.  
  1601.  
  1602.  
  1603.  
  1604. unsigned long GCGchecksum( const char *seq, const long seqlen, unsigned long *checktotal)
  1605. /* GCGchecksum */
  1606. {
  1607.   register long  i, check = 0, count = 0;
  1608.     register char *sp = (char*) seq;
  1609.  
  1610.   for (i = 0; i < seqlen; i++) {
  1611.     count++;
  1612.         check += count * to_upper(*sp);
  1613.       sp++;
  1614.     if (count == 57) count = 0;
  1615.     }
  1616.   check %= 10000;
  1617.   *checktotal += check;
  1618.   *checktotal %= 10000;
  1619.   return check;
  1620. }
  1621.  
  1622. /* Table of CRC-32's of all single byte values (made by makecrc.c of ZIP source) */
  1623. const unsigned long crctab[] = {
  1624.   0x00000000L, 0x77073096L, 0xee0e612cL, 0x990951baL, 0x076dc419L,
  1625.   0x706af48fL, 0xe963a535L, 0x9e6495a3L, 0x0edb8832L, 0x79dcb8a4L,
  1626.   0xe0d5e91eL, 0x97d2d988L, 0x09b64c2bL, 0x7eb17cbdL, 0xe7b82d07L,
  1627.   0x90bf1d91L, 0x1db71064L, 0x6ab020f2L, 0xf3b97148L, 0x84be41deL,
  1628.   0x1adad47dL, 0x6ddde4ebL, 0xf4d4b551L, 0x83d385c7L, 0x136c9856L,
  1629.   0x646ba8c0L, 0xfd62f97aL, 0x8a65c9ecL, 0x14015c4fL, 0x63066cd9L,
  1630.   0xfa0f3d63L, 0x8d080df5L, 0x3b6e20c8L, 0x4c69105eL, 0xd56041e4L,
  1631.   0xa2677172L, 0x3c03e4d1L, 0x4b04d447L, 0xd20d85fdL, 0xa50ab56bL,
  1632.   0x35b5a8faL, 0x42b2986cL, 0xdbbbc9d6L, 0xacbcf940L, 0x32d86ce3L,
  1633.   0x45df5c75L, 0xdcd60dcfL, 0xabd13d59L, 0x26d930acL, 0x51de003aL,
  1634.   0xc8d75180L, 0xbfd06116L, 0x21b4f4b5L, 0x56b3c423L, 0xcfba9599L,
  1635.   0xb8bda50fL, 0x2802b89eL, 0x5f058808L, 0xc60cd9b2L, 0xb10be924L,
  1636.   0x2f6f7c87L, 0x58684c11L, 0xc1611dabL, 0xb6662d3dL, 0x76dc4190L,
  1637.   0x01db7106L, 0x98d220bcL, 0xefd5102aL, 0x71b18589L, 0x06b6b51fL,
  1638.   0x9fbfe4a5L, 0xe8b8d433L, 0x7807c9a2L, 0x0f00f934L, 0x9609a88eL,
  1639.   0xe10e9818L, 0x7f6a0dbbL, 0x086d3d2dL, 0x91646c97L, 0xe6635c01L,
  1640.   0x6b6b51f4L, 0x1c6c6162L, 0x856530d8L, 0xf262004eL, 0x6c0695edL,
  1641.   0x1b01a57bL, 0x8208f4c1L, 0xf50fc457L, 0x65b0d9c6L, 0x12b7e950L,
  1642.   0x8bbeb8eaL, 0xfcb9887cL, 0x62dd1ddfL, 0x15da2d49L, 0x8cd37cf3L,
  1643.   0xfbd44c65L, 0x4db26158L, 0x3ab551ceL, 0xa3bc0074L, 0xd4bb30e2L,
  1644.   0x4adfa541L, 0x3dd895d7L, 0xa4d1c46dL, 0xd3d6f4fbL, 0x4369e96aL,
  1645.   0x346ed9fcL, 0xad678846L, 0xda60b8d0L, 0x44042d73L, 0x33031de5L,
  1646.   0xaa0a4c5fL, 0xdd0d7cc9L, 0x5005713cL, 0x270241aaL, 0xbe0b1010L,
  1647.   0xc90c2086L, 0x5768b525L, 0x206f85b3L, 0xb966d409L, 0xce61e49fL,
  1648.   0x5edef90eL, 0x29d9c998L, 0xb0d09822L, 0xc7d7a8b4L, 0x59b33d17L,
  1649.   0x2eb40d81L, 0xb7bd5c3bL, 0xc0ba6cadL, 0xedb88320L, 0x9abfb3b6L,
  1650.   0x03b6e20cL, 0x74b1d29aL, 0xead54739L, 0x9dd277afL, 0x04db2615L,
  1651.   0x73dc1683L, 0xe3630b12L, 0x94643b84L, 0x0d6d6a3eL, 0x7a6a5aa8L,
  1652.   0xe40ecf0bL, 0x9309ff9dL, 0x0a00ae27L, 0x7d079eb1L, 0xf00f9344L,
  1653.   0x8708a3d2L, 0x1e01f268L, 0x6906c2feL, 0xf762575dL, 0x806567cbL,
  1654.   0x196c3671L, 0x6e6b06e7L, 0xfed41b76L, 0x89d32be0L, 0x10da7a5aL,
  1655.   0x67dd4accL, 0xf9b9df6fL, 0x8ebeeff9L, 0x17b7be43L, 0x60b08ed5L,
  1656.   0xd6d6a3e8L, 0xa1d1937eL, 0x38d8c2c4L, 0x4fdff252L, 0xd1bb67f1L,
  1657.   0xa6bc5767L, 0x3fb506ddL, 0x48b2364bL, 0xd80d2bdaL, 0xaf0a1b4cL,
  1658.   0x36034af6L, 0x41047a60L, 0xdf60efc3L, 0xa867df55L, 0x316e8eefL,
  1659.   0x4669be79L, 0xcb61b38cL, 0xbc66831aL, 0x256fd2a0L, 0x5268e236L,
  1660.   0xcc0c7795L, 0xbb0b4703L, 0x220216b9L, 0x5505262fL, 0xc5ba3bbeL,
  1661.   0xb2bd0b28L, 0x2bb45a92L, 0x5cb36a04L, 0xc2d7ffa7L, 0xb5d0cf31L,
  1662.   0x2cd99e8bL, 0x5bdeae1dL, 0x9b64c2b0L, 0xec63f226L, 0x756aa39cL,
  1663.   0x026d930aL, 0x9c0906a9L, 0xeb0e363fL, 0x72076785L, 0x05005713L,
  1664.   0x95bf4a82L, 0xe2b87a14L, 0x7bb12baeL, 0x0cb61b38L, 0x92d28e9bL,
  1665.   0xe5d5be0dL, 0x7cdcefb7L, 0x0bdbdf21L, 0x86d3d2d4L, 0xf1d4e242L,
  1666.   0x68ddb3f8L, 0x1fda836eL, 0x81be16cdL, 0xf6b9265bL, 0x6fb077e1L,
  1667.   0x18b74777L, 0x88085ae6L, 0xff0f6a70L, 0x66063bcaL, 0x11010b5cL,
  1668.   0x8f659effL, 0xf862ae69L, 0x616bffd3L, 0x166ccf45L, 0xa00ae278L,
  1669.   0xd70dd2eeL, 0x4e048354L, 0x3903b3c2L, 0xa7672661L, 0xd06016f7L,
  1670.   0x4969474dL, 0x3e6e77dbL, 0xaed16a4aL, 0xd9d65adcL, 0x40df0b66L,
  1671.   0x37d83bf0L, 0xa9bcae53L, 0xdebb9ec5L, 0x47b2cf7fL, 0x30b5ffe9L,
  1672.   0xbdbdf21cL, 0xcabac28aL, 0x53b39330L, 0x24b4a3a6L, 0xbad03605L,
  1673.   0xcdd70693L, 0x54de5729L, 0x23d967bfL, 0xb3667a2eL, 0xc4614ab8L,
  1674.   0x5d681b02L, 0x2a6f2b94L, 0xb40bbe37L, 0xc30c8ea1L, 0x5a05df1bL,
  1675.   0x2d02ef8dL
  1676. };
  1677.  
  1678. unsigned long CRC32checksum(const char *seq, const long seqlen, unsigned long *checktotal)
  1679. /*CRC32checksum: modified from CRC-32 algorithm found in ZIP compression source */
  1680. {
  1681.   register unsigned long c = 0xffffffffL;
  1682.   register long n = seqlen;
  1683.     register char *sp = (char*) seq;
  1684.  
  1685.     while (n--) {
  1686.         c = crctab[((int)c ^ (to_upper(*sp))) & 0xff] ^ (c >> 8);
  1687.         sp++;
  1688.         }
  1689.   c= c ^ 0xffffffffL;
  1690.   *checktotal += c;
  1691.   return c;
  1692. }
  1693.  
  1694.  
  1695.  
  1696.  
  1697. short getseqtype( const char *seq, const long seqlen)
  1698. { /* return sequence kind: kDNA, kRNA, kProtein, kOtherSeq, ??? */
  1699.   char  c;
  1700.     register char *sp = (char*) seq;
  1701.   short i, maxtest;
  1702.   short na = 0, aa = 0, po = 0, nt = 0, nu = 0, ns = 0, no = 0;
  1703.  
  1704.   maxtest = (short) min(300, seqlen);
  1705.   for (i = 0; i < maxtest; i++) {
  1706.         c = to_upper(*sp);  sp++;
  1707.     if (strchr(protonly, c)) po++;
  1708.     else if (strchr(primenuc,c)) {
  1709.       na++;
  1710.       if (c == 'T') nt++;
  1711.       else if (c == 'U') nu++;
  1712.       }
  1713.     else if (strchr(aminos,c)) aa++;
  1714.     else if (strchr(seqsymbols,c)) ns++;
  1715.     else if (isalpha(c)) no++;
  1716.     }
  1717.  
  1718.   if ((no > 0) || (po+aa+na == 0)) return kOtherSeq;
  1719.   /* ?? test for probability of kOtherSeq ?, e.g.,
  1720.   else if (po+aa+na / maxtest < 0.70) return kOtherSeq;
  1721.   */
  1722.   else if (po > 0) return kAmino;
  1723.   else if (aa == 0) {
  1724.     if (nu > nt) return kRNA;
  1725.     else return kDNA;
  1726.     }
  1727.   else if (na > aa) return kNucleic;
  1728.   else return kAmino;
  1729. } /* getseqtype */
  1730.  
  1731.  
  1732. void GetSeqStats( const char* seq, const long seqlen, 
  1733.                             short* seqtype, unsigned long* checksum, long* basecount)
  1734. {
  1735.     register unsigned long chk = 0xffffffffL;
  1736.   register long n = seqlen;
  1737.     register char c, *seqp;
  1738.     short         na = 0, aa= 0, no= 0, ns= 0, nt= 0, po= 0, nu= 0;
  1739.     short         seqt = kOtherSeq;
  1740.     unsigned long checks= 0;
  1741.     long     basec= 0;
  1742.     
  1743.     if (seq) {
  1744.         na= aa= no= ns= nt= po= nu= 0;
  1745.         seqp= (char*) seq;
  1746.         while (n--) {
  1747.             c= to_upper(*seqp);
  1748.             seqp++;
  1749.             if (c > ' '  && (c < '0' || c > '9')) { 
  1750.                 chk = crctab[((int)chk ^ c) & 0xff] ^ (chk >> 8);
  1751.             
  1752.                 basec++;
  1753.                 if (strchr( protonly, c)) po++;
  1754.                 else if (strchr(primenuc, c)) {
  1755.                     na++;
  1756.                     if (c == 'T') nt++;
  1757.                     else if (c == 'U') nu++;
  1758.                     }
  1759.                 else if (strchr(aminos, c))  aa++;  
  1760.                 else if (strchr(seqsymbols, c)) ns++;
  1761.                 else if (isalpha(c)) no++;
  1762.                 }
  1763.             }
  1764.       /* checks= chk ^ 0xffffffffL; */
  1765.         if (no > 0 || po+aa+na == 0) seqt = kOtherSeq;
  1766.         else if (po > 0) seqt = kAmino;
  1767.         else if (aa == 0) {
  1768.             if (nu > nt) seqt = kRNA;
  1769.             else seqt = kDNA;
  1770.             }
  1771.         else if (na > aa) seqt = kNucleic;
  1772.         else seqt = kAmino;
  1773.         }
  1774.     if (seqtype) *seqtype= seqt;
  1775.     if (checksum) *checksum= chk ^ 0xffffffffL;
  1776.     if (basecount) *basecount= basec;
  1777. }
  1778.  
  1779.  
  1780. char* compressSeq( const char gapc, const char *seq, const long seqlen, long *newlen)
  1781. {
  1782.   register char *a, *b;
  1783.   register long i;
  1784.   char  *newseq;
  1785.  
  1786.   *newlen= 0;
  1787.   if (!seq) return NULL;
  1788.   newseq = (char*) MALLOC(seqlen+1);
  1789.   if (!newseq) return NULL;
  1790.   for (a= (char*)seq, b=newseq, i=0; *a!=0; a++)
  1791.     if (*a != gapc) {
  1792.       *b++= *a;
  1793.       i++;
  1794.       }
  1795.   *b= '\0';
  1796. #ifdef NOTmacintosh
  1797.     SetPtrSize( (Ptr)newseq, i+1);
  1798. #else
  1799.   newseq = (char*) realloc(newseq, i+1);
  1800. #endif
  1801.   *newlen= i;
  1802.   return newseq;
  1803. }
  1804.  
  1805.  
  1806.  
  1807. /***
  1808. char *rtfhead = "{\\rtf1\\defformat\\mac\\deff2 \
  1809. {\\fonttbl\
  1810.   {\\f1\\fmodern Courier;}{\\f2\\fmodern Monaco;}\
  1811.   {\\f3\\fswiss Helvetica;}{\\f4\\fswiss Geneva;}\
  1812.   {\\f5\\froman Times;}{\\f6\\froman Palatino;}\
  1813.   {\\f7\\froman New Century Schlbk;}{\\f8\\ftech Symbol;}}\
  1814. {\\stylesheet\
  1815.   {\\s1 \\f5\\fs20 \\sbasedon0\\snext1 name;}\
  1816.   {\\s2 \\f3\\fs20 \\sbasedon0\\snext2 num;}\
  1817.   {\\s3 \\f1\\f21 \\sbasedon0\\snext3 seq;}}";
  1818.  
  1819. char *rtftail = "}";
  1820. ****/
  1821.  
  1822.  
  1823.  
  1824.  
  1825.  
  1826.  
  1827. short writeSeq(FILE *outf, const char *seq, const long seqlen,
  1828.                 const short outform, const char *seqid)
  1829. {
  1830.     gFile = outf; gLinestart = 0;
  1831.     return writeSeqCall( readWriteFromFile, seq, seqlen, outform, seqid);
  1832. }
  1833.  
  1834.  
  1835. short writeSeqCall(
  1836.                                 ReadWriteProc writer, 
  1837.                                 const char *seq, const long seqlen,
  1838.                 const short outform, const char *seqid)
  1839. /* dump sequence to standard output */
  1840. {
  1841.   const short kSpaceAll = -9;
  1842. #define kMaxseqwidth  250
  1843.  
  1844.   boolean baseonlynum= false; /* nocountsymbols -- only count true bases, not "-" */
  1845.   short  numline = 0; /* only true if we are writing seq number line (for interleave) */
  1846.   boolean numright = false, numleft = false;
  1847.   boolean nameright = false, nameleft = false, dumb = true;
  1848.   short   namewidth = 8, numwidth = 8;
  1849.   short   spacer = 0, width  = 50, tab = 0;
  1850.   /* new parameters: width, spacer, those above... */
  1851.  
  1852.   short linesout = 0, seqtype = kNucleic;
  1853.   long  i, j, l, l1, ibase;
  1854.   char  idword[31], endstr[30];
  1855.   char  seqnamestore[128], *seqname = seqnamestore;
  1856.   char  s[kMaxseqwidth], *cp;
  1857.   char  nameform[30], numform[30], nocountsymbols[30];
  1858.   unsigned long checksum = 0, checktotal = 0;
  1859.     char    outbuf[512];
  1860.  
  1861. #define FPUTC(c)    {outbuf[0]=c; outbuf[1]=0; (void)(*writer)(outbuf, 255, kRSFile_Write);}
  1862. #define FPUTS(s)    (void)(*writer)(outbuf, 255L, kRSFile_Write)
  1863.  
  1864.     /*fprintf(stderr,"writeSeqCall &writer= %ld  &gFile= %ld\n", writer, gFile);*/
  1865.  
  1866.   gPretty.atseq++;
  1867.   /* skipwhitespace(seqid); //???? */
  1868.   strncpy( seqnamestore, seqid, sizeof(seqnamestore));
  1869.   seqname[sizeof(seqnamestore)-1] = 0;
  1870.  
  1871.   sscanf( seqname, "%30s", idword);
  1872.   sprintf(numform, "%d", seqlen);
  1873.   numwidth= strlen(numform)+1;
  1874.   nameform[0]= '\0';
  1875.  
  1876.   if (strstr(seqname,"checksum") != NULL) {
  1877.     cp = strstr(seqname,"bases");
  1878.     if (cp!=NULL) {
  1879.       for ( ; (cp!=seqname) && (*cp!=','); cp--) ;
  1880.       if (cp!=seqname) *cp=0;
  1881.       }
  1882.     }
  1883.  
  1884.   strcpy( endstr,"");
  1885.   l1 = 0;
  1886.  
  1887.   if (outform == kGCG || outform == kMSF) {
  1888.       /* ! translate quasi-standard '-' indel to '.' that gcg uses */
  1889.       /* must do before checksum... */
  1890.     checksum = GCGchecksum(seq, seqlen, &checktotal);
  1891.     }
  1892.   else
  1893.     checksum = seqchecksum(seq, seqlen, &checktotal);
  1894.  
  1895.   switch (outform) {
  1896.  
  1897.     case kPlain:
  1898.     case kUnknown:    /* no header, just sequence */
  1899.       strcpy(endstr,"\n"); /* end w/ extra blank line */
  1900.       break;
  1901.  
  1902.     case kOlsen:  /* Olsen seq. editor takes plain nucs OR Genbank  */
  1903.     case kGenBank:
  1904.       sprintf(outbuf,"LOCUS       %s       %d bp\n", idword, seqlen);  FPUTS(outbuf);
  1905.       sprintf(outbuf,"DEFINITION  %s  %d bases  %X checksum\n",  
  1906.                         seqname, seqlen, checksum); FPUTS(outbuf);
  1907.    /* sprintf(outbuf,"ACCESSION   %s\n", accnum); FPUTS(outbuf); */
  1908.       sprintf(outbuf,"ORIGIN      \n" ); FPUTS(outbuf);
  1909.       spacer = 11;
  1910.       numleft = true;
  1911.       numwidth = 8;  /* dgg. 1Feb93, patch for GDE fail to read short numwidth */
  1912.       strcpy(endstr, "\n//");
  1913.       linesout += 4;
  1914.       break;
  1915.  
  1916.     case kPIR:
  1917.       /* somewhat like genbank... \\\*/
  1918.       /* sprintf(outbuf,"\\\\\\\n"); FPUTS(outbuf); << only at top of file, not each entry... */
  1919.       sprintf(outbuf,"ENTRY           %s \n", idword); FPUTS(outbuf);
  1920.       sprintf(outbuf,"TITLE           %s  %d bases  %X checksum\n", 
  1921.                     seqname, seqlen, checksum); FPUTS(outbuf);
  1922.    /* sprintf(outbuf,"ACCESSION       %s\n", accnum); FPUTS(outbuf); */
  1923.       sprintf(outbuf,"SEQUENCE        \n" ); FPUTS(outbuf);
  1924.       numwidth = 7;
  1925.       width= 30;
  1926.       spacer = kSpaceAll;
  1927.       numleft = true;
  1928.       strcpy(endstr, "\n///");
  1929.       /* run a top number line for PIR */
  1930.       for (j=0; j<numwidth; j++) FPUTC(' ');
  1931.       for (j= 5; j<=width; j += 5) { sprintf(outbuf,"%10d", j ); FPUTS(outbuf); }
  1932.       FPUTC('\n');
  1933.       linesout += 5;
  1934.       break;
  1935.  
  1936.     case kNBRF:
  1937.       if (getseqtype(seq, seqlen) == kAmino)
  1938.         { sprintf(outbuf,">P1;%s\n", idword );  FPUTS(outbuf);}
  1939.       else
  1940.         { sprintf(outbuf,">DL;%s\n", idword );  FPUTS(outbuf);}
  1941.       sprintf(outbuf,"%s  %d bases  %X checksum\n", 
  1942.                 seqname, seqlen, checksum); FPUTS(outbuf);
  1943.       spacer = 11;
  1944.       strcpy(endstr,"*\n");
  1945.       linesout += 3;
  1946.       break;
  1947.  
  1948.     case kEMBL:
  1949.       sprintf(outbuf,"ID   %s\n", idword ); FPUTS(outbuf);
  1950.   /*  sprintf(outbuf,"AC   %s\n", accnum ); FPUTS(outbuf); */
  1951.       sprintf(outbuf,"DE   %s  %d bases %X checksum\n", 
  1952.                         seqname, seqlen, checksum); FPUTS(outbuf);
  1953.       sprintf(outbuf,"SQ             %d BP\n", seqlen ); FPUTS(outbuf);
  1954.       strcpy(endstr, "\n//"); /* 11Oct90: bug fix*/
  1955.       tab = 4;     /** added 31jan91 */
  1956.       spacer = 11; /** added 31jan91 */
  1957.       width = 60;
  1958.       linesout += 4;
  1959.       break;
  1960.  
  1961.     case kGCG:
  1962.       sprintf(outbuf,"%s\n", seqname ); FPUTS(outbuf);
  1963.    /* sprintf(outbuf,"ACCESSION   %s\n", accnum );  FPUTS(outbuf);*/
  1964.       sprintf(outbuf,"    %s  Length: %d  (today)  Check: %d  ..\n", 
  1965.                 idword, seqlen, checksum ); FPUTS(outbuf);
  1966.       spacer = 11;
  1967.       numleft = true;
  1968.       strcpy(endstr, "\n");  /* this is insurance to help prevent misreads at eof */
  1969.       linesout += 3;
  1970.       break;
  1971.  
  1972.     case kStrider: /* ?? map ?*/
  1973.       sprintf(outbuf,"; ### from DNA Strider ;-)\n" ); FPUTS(outbuf);
  1974.       sprintf(outbuf,"; DNA sequence  %s  %d bases  %X checksum\n;\n", 
  1975.                 seqname, seqlen, checksum); FPUTS(outbuf);
  1976.       strcpy(endstr, "\n//");
  1977.       linesout += 3;
  1978.       break;
  1979.  
  1980.     case kFitch:
  1981.       sprintf(outbuf,"%s  %d bases  %X checksum\n", 
  1982.                 seqname, seqlen, checksum); FPUTS(outbuf);
  1983.       spacer = 4;
  1984.       width = 60;
  1985.       linesout += 1;
  1986.       break;
  1987.  
  1988.     case kPhylip2:
  1989.     case kPhylip4:
  1990.       /* this is version 3.2/3.4 -- simplest way to write
  1991.         version 3.3 is to write as version 3.2, then
  1992.         re-read file and interleave the species lines */
  1993.       if (strlen(idword)>10) idword[10] = 0;
  1994.       sprintf(outbuf,"%-10s  ", idword ); FPUTS(outbuf);
  1995.       l1  = -1;
  1996.       tab = 12;
  1997.       spacer = 11;
  1998.       break;
  1999.  
  2000.     case kASN1:
  2001.       seqtype= getseqtype(seq, seqlen);
  2002.       switch (seqtype) {
  2003.         case kDNA     : cp= "dna"; break;
  2004.         case kRNA     : cp= "rna"; break;
  2005.         case kNucleic : cp= "na"; break;
  2006.         case kAmino   : cp= "aa"; break;
  2007.         case kOtherSeq: cp= "not-set"; break;
  2008.         }
  2009.       sprintf(outbuf,"  seq {\n" ); FPUTS(outbuf);
  2010.       sprintf(outbuf,"    id { local id %d },\n", gPretty.atseq ); FPUTS(outbuf);
  2011.       sprintf(outbuf,"    descr { title \"%s\" },\n",  seqid ); FPUTS(outbuf);
  2012.       sprintf(outbuf,"    inst {\n", dumb ); FPUTS(outbuf);
  2013.       sprintf(outbuf,"      repr raw, mol %s, length %d, topology linear,\n", 
  2014.                 cp, seqlen ); FPUTS(outbuf);
  2015.       sprintf(outbuf,"      seq-data\n" ); FPUTS(outbuf);
  2016.       if (seqtype == kAmino)
  2017.         { sprintf(outbuf,"        iupacaa \""); FPUTS(outbuf);}
  2018.       else
  2019.         { sprintf(outbuf,"        iupacna \"" ); FPUTS(outbuf);}
  2020.       l1  = 17;
  2021.       spacer = 0;
  2022.       width  = 78;
  2023.       tab  = 0;
  2024.       strcpy(endstr,"\"\n      } } ,");
  2025.       linesout += 7;
  2026.       break;
  2027.  
  2028.     case kPAUP:
  2029.       nameleft= true;
  2030.       namewidth = 9;
  2031.       spacer = 21;
  2032.       width  = 100;
  2033.       tab  = 0; /* 1; */
  2034.       /* strcpy(endstr,";\nend;"); << this is end of all seqs.. */
  2035.       /* do a header comment line for paup */
  2036.       sprintf(outbuf,"[Name: %-16s  Len:%6d  Check: %8X]\n", 
  2037.                 idword, seqlen, checksum ); FPUTS(outbuf);
  2038.       linesout += 1;
  2039.       break;
  2040.  
  2041.     case kPretty:
  2042.       numline= gPretty.numline;
  2043.       baseonlynum= gPretty.baseonlynum;
  2044.       namewidth = gPretty.namewidth;
  2045.       numright = gPretty.numright;
  2046.       numleft = gPretty.numleft;
  2047.       nameright = gPretty.nameright;
  2048.       nameleft = gPretty.nameleft;
  2049.       spacer = gPretty.spacer + 1;
  2050.       width  = gPretty.seqwidth;
  2051.       tab  = gPretty.tab;
  2052.       /* also add rtf formatting w/ font, size, style */
  2053.       if (gPretty.nametop) {
  2054.         sprintf(outbuf,"Name: %-16s  Len:%6d  Check: %8X\n", 
  2055.                     idword, seqlen, checksum ); FPUTS(outbuf);
  2056.         linesout++;
  2057.         }
  2058.       break;
  2059.  
  2060.     case kMSF:
  2061.       sprintf(outbuf," Name: %-16s Len:%6d  Check:%5d  Weight:  1.00\n",
  2062.                      idword, seqlen, checksum ); FPUTS(outbuf);
  2063.       linesout++;
  2064.       nameleft= true;
  2065.       namewidth= 15; /* need MAX namewidth here... */
  2066.       sprintf(nameform, "%%+%ds ",namewidth);
  2067.       spacer = 11;
  2068.       width  = 50;
  2069.       tab = 0; /* 1; */
  2070.       break;
  2071.  
  2072.     case kIG:
  2073.       sprintf(outbuf,";%s  %d bases  %X checksum\n",  
  2074.                 seqname, seqlen, checksum ); FPUTS(outbuf);
  2075.       sprintf(outbuf,"%s\n",  idword ); FPUTS(outbuf);
  2076.       strcpy(endstr,"1"); /* == linear dna */
  2077.       linesout += 2;
  2078.       break;
  2079.  
  2080.     default :
  2081.     case kZuker: /* don't attempt Zuker's ftn format */
  2082.     case kPearson:
  2083.       sprintf(outbuf,">%s  %d bases  %X checksum\n",  
  2084.                 seqname, seqlen, checksum ); FPUTS(outbuf);
  2085.       linesout += 1;
  2086.       break;
  2087.     }
  2088.  
  2089.   if (*nameform==0) sprintf(nameform, "%%%d.%ds ",namewidth,namewidth);
  2090.   if (numline) sprintf(numform, "%%%ds ",numwidth);
  2091.   else sprintf(numform, "%%%dd ",numwidth);
  2092.   strcpy( nocountsymbols, kNocountsymbols);
  2093.   if (baseonlynum) {
  2094.     if (strchr(nocountsymbols,gPretty.gapchar)==NULL) {
  2095.       strcat(nocountsymbols," ");
  2096.       nocountsymbols[strlen(nocountsymbols)-1]= gPretty.gapchar;
  2097.       }
  2098.     if (gPretty.domatch && (cp=strchr(nocountsymbols,gPretty.matchchar))!=NULL) {
  2099.       *cp= ' ';
  2100.       }
  2101.     }
  2102.  
  2103.   if (numline) {
  2104.    *idword= 0;
  2105.    }
  2106.  
  2107.   width = min(width,kMaxseqwidth);
  2108.   for (i=0, l=0, ibase = 1; i < seqlen; ) {
  2109.  
  2110.     if (l1 < 0) l1 = 0;
  2111.     else if (l1 == 0) {
  2112.       if (nameleft) { sprintf(outbuf, nameform,  idword ); FPUTS(outbuf);}
  2113.       if (numleft) { if (numline) { sprintf(outbuf, numform,  "" ); FPUTS(outbuf);}
  2114.                     else { sprintf(outbuf, numform, ibase );  FPUTS(outbuf);}
  2115.                                         }
  2116.       for (j=0; j<tab; j++) FPUTC(' ');
  2117.       }
  2118.  
  2119.     l1++;                 /* don't count spaces for width*/
  2120.     if (numline) {
  2121.       if (spacer==kSpaceAll || (spacer != 0 && (l+1) % spacer == 1)) {
  2122.         if (numline==1) FPUTC(' ');
  2123.         s[l++] = ' ';
  2124.         }
  2125.       if (l1 % 10 == 1 || l1 == width) {
  2126.         if (numline==1) { sprintf(outbuf,"%-9d ", i+1 ); FPUTS(outbuf);}
  2127.         s[l++]= '|'; /* == put a number here */
  2128.         }
  2129.       else s[l++]= ' ';
  2130.       i++;
  2131.       }
  2132.  
  2133.     else {
  2134.       if (spacer==kSpaceAll || (spacer != 0 && (l+1) % spacer == 1))
  2135.         s[l++] = ' ';
  2136.       if (!baseonlynum) ibase++;
  2137.       else if (0==strchr(nocountsymbols,seq[i])) ibase++;
  2138.       s[l++] = seq[i++];
  2139.       }
  2140.  
  2141.     if (l1 == width || i == seqlen) {
  2142.       if (outform==kPretty) for ( ; l1<width; l1++) {
  2143.         if (spacer==kSpaceAll || (spacer != 0 && (l+1) % spacer == 1))
  2144.           s[l++] = ' ';
  2145.         s[l++]=' '; /* pad w/ blanks */
  2146.         }
  2147.       s[l] = '\0';
  2148.       l = 0; l1 = 0;
  2149.  
  2150.       if (numline) {
  2151.         if (numline==2) { sprintf(outbuf,"%s", s );  FPUTS(outbuf);} /* finish numberline ! and | */
  2152.         }
  2153.       else {
  2154.         if (i == seqlen) { sprintf(outbuf,"%s%s", s,endstr ); FPUTS(outbuf); }
  2155.         else { sprintf(outbuf,"%s", s ); FPUTS(outbuf);}
  2156.         if (numright || nameright) FPUTC(' ');
  2157.         if (numright)  { sprintf(outbuf,numform,  ibase-1 ); FPUTS(outbuf);}
  2158.         if (nameright) { sprintf(outbuf, nameform, idword ); FPUTS(outbuf);}
  2159.         }
  2160.       FPUTC('\n');
  2161.       linesout++;
  2162.       }
  2163.     }
  2164.   return linesout;
  2165.     
  2166. #undef FPRINTF 
  2167. #undef FPUTC 
  2168. }  /*writeSeq*/
  2169.  
  2170.  
  2171.  
  2172. /* End file: ureadseq.c */
  2173.